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Abstract 

We present the results from a series of lattice simulations of the 
charmonium system using a highly-improved NRQCD action, both in 
the quenched approximation, and with nt = 2 light dynamical quarks. 
The spectra show some evidence for quenching effects of roughly 10% 
in the S- and P-hyperfine spin splittings — probably too small to ac- 
count for the severe underestimates in these quantities seen in previ- 
ous quenched charmonium simulations. We also find estimates for the 
magnitude of other systematic effects — in particular, the choice of the 
tadpole factor can alter spin splittings at the 10-20% level, and 0(a s ) 
radiative corrections may be as large as 40% for charmonium. Wc 
conclude that quenching is just one of a collection of important effects 
that require attention in precision heavy-quark simulations. 

1 Charmonium on the lattice 

One of the most rapidly expanding sectors of lattice QCD in the last decade 
has been the study of heavy-quark systems. Lattice simulations have suc- 
cessfully reproduced the broad structure of the heavy hadron spectrum, 
providing a solid piece of evidence for the correctness of QCD. Discrepan- 
cies at the level of the hyperfine structure still persist however, and in some 
cases these are uncomfortably large. 

This paper describes a series of highly-improved non-relativistic simula- 
tions of the charmonium system, with the aim of estimating the sizes of var- 
ious systematic uncertainties influencing the spectrum. An understanding 
of the relative influence of these uncertainties on the heavy-quark spectrum 
is vital to the interpretation of the current state of lattice simulations. 

One very successful approach to simulating heavy quark systems utilises 
the NRQCD formalism |2), where the quark dynamics are governed by 
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an effective non-reiativistic Hamiltonian, expanded in powers of the heavy- 
quark velocity. For the bottom and charm quarks, v 2 ~ 0.1 and v 2 ~ 0.3 
respectively, and so we expect to achieve some success with a non-relativistic 
theory. Simulations of heavy-light and heavy-heavy charm and bottom sys- 
tems have shown that NRQCD captures much of the correct physics of the 
heavy quarks. Understanding the remaining systematic errors in heavy- 
quark simulations has become a major focus of the lattice NRQCD commu- 
nity. 

The first report of a high-statistics NRQCD simulation of charmonium 
appeared in 1995 by Davies et al. Q. The authors used a NRQCD Hamilto- 
nian with relativistic and discretisation errors corrected to 0(v 4 ) to measure 
ground and excited S, P and D state energies in the quenched approxima- 
tion. Agreement with experiment was very promising, with discrepancies at 
the order of 10-30% in S- and P-state hyperfine spin-splittings, in agreement 
with the expected size of the next-order corrections. 

Disturbingly, charmonium simulations incorporating 0{v & ) corrections 
[H] showed a large decrease in hyperfine spin splittings, taking theoretical 
predictions as much as 50% further away from experimental values. These 
simulations also demonstrated a large dependence on the definition of the 
tadpole correction factor. Given the size of v 2 for charmonium, sizeable 
0(v e ) corrections are not surprising; however, the disappointingly large dis- 
crepancies in the spectrum with such a highly improved theory give pause 
to the future of charmonium simulations. Evidently, the NRQCD expansion 
converges slowly for the charm quark. 

Even in the less-relativistic T system, the same highly-improved NRQCD 
action has not provided conclusive agreement with experiment [||, g. Cer- 
tainly, NRQCD to C(f 6 ) is not a closed problem. 

The difficulties with the hyperfine spectrum are not limited to the NRQCD 
approach. A report on the status of charmonium simulations with the rel- 
ativistic Fermilab approach in 1993 || cited a 20-30% shortfall for the S- 
state hyperfine splitting using an SW-improved Wilson action. In 1999, 
the UKQCD collaboration reported on a tadpole- and SW-improved sim- 
ulation of charmonium ||; their results for the S'-hyperfine splitting were 
roughly 40% below the experimental value. Both of these simulations used 
the quenched approximation, and the inclusion of dynamical quark loops 
would increase the hyperfine splittings. In the 1993 report, quenching ef- 
fects were estimated to be as large as 40%, however this seems optimistic — 
corrections at the 5-15% level seem more typical in full QCD simulations of 
both the T system [|j and of light hadrons [10]. 

A very recent report from the CP-PACS collaboration |?j describes un- 
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quenched simulations of charmonium and bottomonium using NRQCD over 
a range of lattice spacings and sea-quark masses, with rif = 2 SW-improved 
light sea-quarks. In that report, the authors concentrate mostly on sim- 
ulations of the 66-system, though some charmonium results are presented. 
Their results indicate a significant increase in the S-state hyperfine split- 
tings as the sea-quark mass decreases towards the chiral limit, though no 
effect is seen on the P-states. 

We have performed a series of highly-improved NRQCD simulations to 
examine the various systematic uncertainties influencing the charmonium 
spectrum. We first concentrate on the effects of dynamical quark loops. 
If these account for the majority of the hyperfine splitting discrepancy in 
charmonium then we expect to find a large increase in the splittings when 
dynamical quarks are included, even in the NRQCD formalism. We examine 
this effect using an ensemble of unquenched configurations provided by the 
MILC collaboration, seeking to establish whether the effects of dynamical 
quarks are sufficient to reconcile the hyperfine discrepancy. This work does 
not aim to provide the definitive unquenched charmonium spectrum. 

The remainder of the paper is devoted to an examination of other sys- 
tematic effects. Simulations with two common definitions of the tadpole 
correction factor result in significantly different spectra, and we find a rough 
estimate of the effect of 0{a s ) radiative corrections to the NRQCD expan- 
sion coefficients. Finally, we note a sizeable shift in the hyperfine splittings 
due to an instability in the standard form for the heavy-quark propagator's 
evolution equation. Each of these effects is contrasted with the estimated 
magnitude of the unquenching error, which leads us to several conclusions 
about NRQCD simulations of charmonium, and heavy-quark simulations in 
general. 



2 The standard lattice NRQCD formalism 

The NRQCD Hamiltonian is typically presented as an expansion in powers 
of the heavy-quark velocity. A highly-improved NRQCD Hamiltonian, with 
corrections to 0(v 6 ) in the velocity expansion Q, is 

H = H + 5H V 4 + 5H v e , (I) 

where 

-AW 
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is the leading kinetic Schrodinger operator, and the 0(v A ) and C(t; 6 ) cor- 
rections are 

2 
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5H v6 = -c 7 ^{A^,a.B} 
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{A( 2 ),cj- (AxE-ExA)} 



A iiide signifies the use of improved versions of the lattice operators that 
remove the leading discretisation errors: the improved lattice derivatives A 
and A^ 2 ) are given by 

a 2 

A^(n) = A^(n)-yA^(n), 

A^(n) = A 2 ^n) + ^(A 2 )^(n), (5) 

while the fields E{ = F44 and Bi = \(-ijkFjk are taken from an improved 
gauge field tensor Q |2j , 

F„u(n) = ~i^(n) - - [u^{n)F^(n + £)E/J(n) + C/J(n - A)-^(" " " A) 

-^(n)F^(n + />)C/t( n ) + [/t(n - *>)i^,(n - u)U u {n - v)] . (6) 

All lattice operators are tadpole improved [O], by dividing all instances of 
the link operators U^(n) by the tadpole correction factor no, 

no 

This means, for example, that the gauge E and B fields are adjusted by a 
factor of Uq. Much evidence exists for the superiority of the Landau defini- 
tion of the tadpole factor, 

,l n. 



u^ = {-TtU,) , (8) 
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over the plaquette definition, 

/ 1 \ l / 4 
4 = {-TrP^ . (9) 

For example, Uq leads to smaller corrections to hyperfine splittings, and 
better scaling of quarkonium masses 0, 0|; it restores rotational invariance 
to a greater degree in the static quark potential [Oj; and it results in closer 
agreement between the tadpole-improved value and the perturbative value 
for the 'clover' coefficient c sw in the Sheikholeslami-Wohlert action [111]. We 
have used both the Landau and plaquette definition in our simulations. 

Since the quarks and antiquarks are decoupled in the non-relativistic 
theory, the heavy-quark Green's function may be found from an evolution 
equation, 

G t+1 = (1-^^(1-^(1- aSH) G t , (10) 
with the initial time-step given by 

The (1 — aH) factors are linear approximations to the continuum evolution 
operator e . The 'stabilisation parameter' s appearing in Equations (|l|) 
and ( JlOj ) improves the approximation to the time evolution operator e aH . 

To complement the use of a highly-improved quark Hamiltonian, we use 
a tadpole and 'rectangle' improved action for the gauge fields [O, 

s G = -pJ2 ( A^VM - + R^)) > (i 2 ) 

where P MJ/ (n) and -R^^ represent the traces of 1 x 1 plaquettes and 2x1 
rectangles of link operators respectively. 

Operators for the various quarkonium states have the form 

M(t) = J2^(n,t)T(n)xHn,t), (13) 

n 

where ijr and are t ne quark and antiquark creation operators, and T(n) 
provides the appropriate spin and spatial wavefunction quantum numbers. 
Operators for the lowest-lying S, P and D states are given in a number 
of references ||, f|; using these, we have constructed propagators for each 
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of the 2S+1 Lj = 1 S , 3 5i, 1 Pi, 3 P , 3 -Pi and 3 P 2 states. Only one spin 
polarisation of each of the triplet states was used. 

To reduce the effects of excited-state contamination and improve the 
operators' overlap with the true meson ground-state wavefunctions, we have 
used a gauge-invariant smearing function, replacing 

I» -> T(n)(f> sm {n) . (14) 

A simple and effective choice for (j) sm is ]T5| 

sm (e,n s ) = (l + eA 2 )™\ (15) 

The weighting factor e and number of smearing iterations n s were tuned to 
optimise the overlap with the ground state. 



3 Details of the simulations 

We have performed a number of different simulations of the charm system, to 
compare the magnitudes of various systematic effects on the spectrum. We 
obtained results with the NRQCD Hamiltonian in Equation (jl|) truncated 
to 0(v 4 ) and 0(v 6 ), with both the Landau and plaquette definitions for the 
tadpole factor uq. 

To examine the size of dynamical quark effects, we obtained an ensemble 
of 200 unquenched gauge field configurations, generously provided by the 



MILC collaboration [16]. The configurations were created with the Wilson 
gluon action at /3 = 5.415, with two flavours of staggered dynamical quarks 
at m = 0.025. This light quark mass corresponds to a pseudoscalar-to- 
vector meson mass ratio of m ps /m v ~ 0.45. The lattice volume of these 
configurations is 16 3 x 32 — with a spacing of a ~ 0.16 fm (determined from 
the charmonium spectrum as described below), this corresponds to a lattice 
extending roughly 2.5 fermi in each spatial direction. 

We produced an ensemble of quenched configurations with both the Lan- 
dau and plaquette tadpole definitions, using the improved action in Equa- 
tion [l^. We found that, using Landau and plaquette tadpoles respectively, 
(3 = 2.1 and (3 = 2.52 give almost the same lattice spacing as the unquenched 
configurations. These results agree with the spacings given in Reference @] 
at the same values of (3. We created 100 configurations in each case, with 
lattice volume 12 3 x 24, the largest we were able to manage with our com- 
putational resources. Given the small physical size of the heavy mesons, 
however, the difference in volume between the quenched and unquenched 
configurations should not have an effect on our results. 
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The lattice spacing was determined for each ensemble using the spin- 
averaged PS splitting, for charmonium E(P — S) = 458 MeV. This splitting 
is known to be quite independent of the heavy quark mass, falling only 
slightly to 440 MeV for bottomonium, and so serves as a stable quantity for 
determining the physical lattice spacing. We have collected the parameters 
of our simulations together in Table |]. 

The kinetic mass of a boosted state with momentum p is defined by 

E(P) = E(0) + ^f k + 0(V A ) ■ (16) 

The bare charm quark mass Mq is tuned by requiring that the kinetic mass 
of the 1 So charmonium state agrees with the experimental mass of the r] c , 
M Vc = 2.98 GeV. We created correlators for a boosted state with p = 
0, 0), where L is the spatial extent of the lattice. The tuned bare masses, 
and their corresponding physical (kinetic) masses for the 1 Sq, are shown in 
Table |. 

Meson correlators were calculated for the various charmonium states, 
using smeared meson operators with n s = 8 and e = 1/12 in Equation ( |i~5|) at 
both the source and sink. To decrease statistical uncertainties, we calculated 
more than one meson correlator for each gauge field configuration. Meson 
sources were situated at four different spatial origins — (0, 0, 0), (L/2, L/2, 0), 
(L/2,0, L/2) and (0, L/2, L/2) — and starting from two time slices, at t = 
and t = 12, for a total of 800 meson correlator measurements for each state. 

Statistical correlations will exist between the multiple measurements of 
the propagators within each configuration, however the small size of QQ sys- 
tems (the cc is roughly 0.5 fm in radius) is some justification for this practice. 
The correlations are expected to be small, as noted in other charmonium 
studies with similar lattice spacings j|, [| . 

Masses for the various cc states were found by fitting the correlators with 
a single exponential, 

G M {t>t min ) = c M e- EMt (17) 

after a minimum time t m i n , allowing for suitable suppression of excited state 
contributions. Energy splittings between correlated states, such as the S- 
state hyperfme splitting AE = E( s Si) — E( 1 So), can often be extracted 
more precisely by fitting to a ratio of their two correlators, 

Vm - GB{t) -> CB6 ' EBt - C B^ (EA+SE)t _ r -5Et (m) 
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We used ratio fits to extract the 5-state hyperfine splitting, and the kinetic 
mass from the boosted So state. Attempts to extract P-state hyperfine 
splittings in this manner were unsuccessful, as statistical noise overtook the 
very small signal before a reasonable plateau emerged. Single-exponential 
fits, however, resolved the three 3 P levels. We have not employed a bootstrap 
analysis for the fit results, which may suggest we have overestimated the 
statistical uncertainties. 

In the following sections, we present the results for a range of simula- 
tions, incorporating all combinations of quenched and unquenched gauge 
configurations, 0(v 4 ) and 0(v e ) correction terms, and Landau and plaque- 
tte tadpole factors. 

3.1 Quenched results 

An example of the quality of the correlator data is shown in Figures I and 
||, plots of the effective masses for the 1 5'o, Pi and 3 Po from the simulation 
using the Landau tadpole factor. The meson propagators were fit with single 
exponentials over a range of time intervals (t m i n : t max ). An indication of 
the convergence of these fits is given in Table ||, where the fit results are 
shown for the 0(v e ) simulations using the plaquette tadpole factor. The 
results presented in this table are representative of all of the charmonium 
spectra we present here. The two S-states had a much cleaner signal than 
the four P-states, evident in the lower value for t max used for the P-state 
fits. 

Table ^ contains the final results for the quenched charmonium mass 
fits. We considered the ground-state for each meson propagator to have 
properly emerged when three consecutive t m i n : t max intervals gave results 
that agreed within statistical errors; the meson mass was then taken as the 
middle of these three values. The masses are given in both lattice units and 
physical units, using the values for a -1 in Table |l] to provide the physical 
energy scale. The simulated spectra are displayed in Figures ^ and ||], shown 
against the experimental data. 

3.2 Unquenched Results 

Given the similar lattice spacings of the MILC configurations and our own 
quenched ensembles, we have used almost the same parameter set for the 
unquenched charmonium simulations — the lower half of Table |] shows the 
specific parameters used. The results of the unquenched simulations are 
contained in Table ||, with the physical energy scale set by a -1 = 1.21(2) 
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GeV, again from the spin-averaged P-S splitting. The spectra are shown in 
Figures || and ^. 

4 Discussion of the Spectra 

A cursory comparison of the quenched and unquenched results shows that, 
while the qualitative structure of the spectrum appears, precision NRQCD 
simulations of the charmonium system have a number of issues yet to be 
resolved. This is most readily seen in the hyperfine splittings, which are 
collected in Figures |7| and || and compared in Table |H[ 

Consider first the quenched results. The 0(v 6 ) corrections lead to a 
disturbingly large decrease in the hyperfine splittings, taking them further 
away from the experimental values by as much as 60%. The situation for 
the plaquette-tadpole simulations is strikingly bad, where the 3 P states 
appear in the wrong order. This reversal is corrected in the Landau-tadpole 
simulations, though the hyperfine splittings are still badly underestimated. 

These difficulties are not new — Trottier j|] first drew attention to the 
large 0(v & ) corrections to the 5-state hyperfine splitting in 1996, and noted 
a possible problem with the 3 P-state ordering. Trottier and Shakespeare 
[ |12| | examined the effects of the different tadpole definitions Uq and Uq on 
the <S-state hyperfine splitting. They performed C(u 6 )-improved NRQCD 
simulations using both tadpole schemes, across a wide range of lattice spac- 
ings, and drew a number of important conclusions; most notably, the 0(v & ) 
hyperfine corrections with Landau tadpoles were significantly smaller than 
the plaquette tadpole results. 

We have confirmed a number of these results here, and in particular 
clearly resolved the extremely poor 3 P-state behaviour, most notably when 
Uq is used. This may simply be a problem due to the bare charm mass 
falling below one in these simulations. However, the Uq simulations lead 
to a higher bare c-quark mass for a given lattice spacing, and the very low 
P-state hyperfine splitting even with aMo > 1 suggests that these problems 
extend beyond the size of the bare mass. 

4.1 Evidence for Quenching Effects? 

The large discrepancies in spin-dependent splittings would be less worrisome 
if quenching were seen to have a considerable effect on the spectrum, as sug- 
gested in ||. Sadly, this does not seem to be the case. There is some evidence 
for a difference between the quenched and unquenched simulations in the 
0(v e ) S'-hyperfine data, perhaps as much as ten percent. However, given 
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the apparent size of other systematic uncertainties, no great significance can 
be attached to these differences. 

We must address the difference between the quenched and unquenched 
gluon actions — the MILC configurations were created with the Wilson pla- 
quette action, while we have employed the rectangle-improved action for the 
quenched lattices. We therefore anticipate an 0{a) error entangled with the 
effects of the dynamical quarks. Our quenched 0(v 4 ) results can be com- 
pared with the results from Reference ||, where the plaquette action was 
used at roughly the same lattice spacing. We see a ~ 10 MeV difference 
between the S'-hyperfine splittings in the two simulations. 

We wish to reiterate our goal, however, to see whether the dynamical 
quark effects are large or small. The S-hyperfine splitting, even in relativistic 
simulations, falls short of experiment by 40 to 50 MeV. An unquenching 
effect of this magnitude would be visible, even taking differences in gluon 
action into account. No such effect was observed in these simulations, and 
we therefore suggest that quenching effects are small in this sense. 

This conclusion is supported by results in high-precision T simulations 
|0, |6|, where the P-state hyperfine splitting is still somewhat underestimated 
in unquenched simulations of this highly-nonrelativistic system, despite the 
use of the 0(v 6 ) improved NRQCD action. Very recently, a 10% sea-quark 
effect was seen in the hyperfine splittings of the charmonium and bottomo- 
nium system in Reference {?J, but differences between the nj = and nj = 2 
P-state splittings were not significant compared with other systematic un- 
certainties. Recent results with unquenched lattices in the B meson spec- 
trum have also shown no significant differences between nf = and nj = 2 
dynamical quark flavours 

4.2 Other Systematic Errors 

The preceding results suggest that agreement between lattice simulations 
and experiment in quarkonium systems will likely not improve through the 
effects of dynamical quarks alone. In the remainder of this section we explore 
various other systematic errors that impact on heavy-quark simulations, as 
a contrast to the small quenching effects found above. 

4.2.1 The Choice of the Tadpole Factor 

We have seen, as others have previously, large differences between results 
using the Landau tadpole factor Uq , and those with the plaquette definition 
Uq. In our own simulations, the size of the 0(v e ) corrections with Uq are 
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significantly smaller than the plaquette tadpole results. This is not sur- 
prising: the E and B fields are each multiplied by a factor of tig 4 in the 
tadpole-improved theory. On our lattices, 

Uo_Y _ j 1-24 (Quenched) , . 

u l) ~ \ 1.30 (MILC) Uyj 

Terms in the NRQCD Hamiltonian linear in E or B will differ by as much 
as 30% between the different tadpole improvement schemes. 

As noted earlier, the evidence in favour of Landau tadpoles is strong. Our 
simulations offer further support, particularly in the 3 P-state behaviour, 
though the more salient issue here is that tadpole effects are at least as 
important as quenching effects in our simulations. 



4.2.2 Radiative Corrections 

We expect some effect on the spectrum from high-momentum modes that 
are cut off by the finite lattice spacing. These high-energy effects may be 
calculated in perturbative QCD as 0(a s ) radiative corrections to the coef- 
ficients of the NRQCD expansion, and there are indications that these may 
be large for the charm quark. Lattice perturbation theory calculations of 
corrections to c\ and C5, the 'kinetic' terms in Equation (|l]), have been com- 
pleted by Morningstar |fTq| . The corrections are roughly 10% or less for the 
bottom quark, but rise dramatically as the bare quark mass falls below one 
(in lattice units). In typical simulations, the bare charm quark mass sits 
close to unity, and so these corrections may become quite significant. 

It is possible to find these radiative corrections without performing long 
calculations in lattice perturbation theory, by using Monte Carlo simulations 
at very high values of []l9| . Such 'non-perturbative' perturbative results 
have been obtained by Trottier and Lepage [^] for the spin-dependent C4 
term in the 0{v A ) NRQCD Hamiltonian, Equation (|l]). Unfortunately, ra- 
diative corrections to the remaining terms in the NRQCD Hamiltonian have 
not been calculated to date. 

We performed a 'toy' simulation to roughly estimate the effects of 0(a s ) 
corrections to all terms in the NRQCD Hamiltonian, replacing the tree level 
coefficients q = 1 with Ci = 1 ± a s . A rough estimate of a s can be made 
from the (tadpole-improved) parameters of our simulations, 

a<w <.)-«ri + o(» 2 )^= i ^. (20) 
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For our values of (5 and uq, this gives a s ~ 0.15-0.2. For the three terms 
in the Hamiltonian where perturbative analysis has been performed, we 
used the calculated values (l^, pQj ; for the remaining terms, we varied the 
coefficients between 0.8 and 1.2. 

Altering the coefficients in this way, we found that the charmonium 
S- and P-hyperfine splittings changed by as much as 10-40%, depending 
on the sign of the corrections for each individual Cj. While this is only 
a crude estimate, it is clear that the effects of radiative corrections may 
be as important as quenching effects for heavy-quark systems. Accurate 
determinations of the remaining 0(a s ) corrections are sorely needed. 

4.2.3 Improving the Evolution Equation 

The evolution equation we presented in Section || for the heavy-quark prop- 
agator, Equation (|H]), contains better-than-linear approximations to the 
exponential e Ht for the terms involving the zeroth-order Hamiltonian Hq, 
but only a linear approximation for the correction terms 5H. Noting that 
the high-order corrections are quite large for charmonium, it is conceivable 
that this lowest-order approximation is too severe. A similar conclusion was 
made by Lewis and Woloshyn of their NRQCD simulations of the D meson 
spectrum pj| . The authors were able to remove some spurious effects due 
to large vacuum expectation values of one of the high-order terms in their 
NRQCD Hamiltonian p3], by improving the exponential approximation for 
the 5H terms in the evolution equation. 

The coefficients of the 0(v 6 ) terms include high powers of Mq 1 and Uq 1 , 
and it is conceivable that for the charm quark, with aMo ~ 1, the (l—a5H v e ) 
approximation is poor. We examined this possibility for the 0(v e ) terms, 
by using an improved form for the evolution equation that incorporates a 
'stabilisation' parameter for the correction terms, with the replacement 

(l-o*jff)->(l-^)'\ (21) 

We have performed a simulation with this alteration to the evolution 
equation, with s$ = 4. Otherwise, all other parameters were kept the same 
as the previous Landau-tadpole quenched simulations. In general, altering 
the evolution equation will lead to a change in the bare charm quark mass 
Mo. In this case we found that Mq = 1.15 once again gave a value of 3.0(1) 
GeV for the 1 S'o mass. 

The the improved evolution equation altered the S-hyperfine splitting 
significantly, increasing it by roughly 40% to 70 MeV. The statistical uncer- 
tainties in the P-hyperfine splittings were large, though a similar increase 
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seems likely. These results suggest the linear approximation (1 — a8H) typ- 
ically used in NRQCD simulations is not sufficiently accurate for the large 
corrections encountered at the charm quark mass. 

5 Conclusions 

Lattice NRQCD simulations of heavy-quark systems have evolved greatly 
over the last decade. By incorporating high-order interaction terms to 
counter relativistic and discretisation errors, simulations now routinely pro- 
duce results that agree with experiment at the 10-30% level. However, 
stubborn discrepancies remain in highly-improved simulations, typically per- 
formed in the quenched approximation, or at tree-level in the 0{a s ) expan- 
sion, or both. To proceed further, all remaining systematic errors must be 
addressed. 

Past studies strongly suggest that the NRQCD expansion converges 
slowly for the charm quark, with the leading and next-to-leading order 
corrections apparently oscillating in sign. To 0(v e ), the hyperfine spin- 
splittings fall short of experimental values by 50% or more. Without know- 
ing the magnitude of the next-order corrections in the velocity expansion, 
the question of reducing the disparities in the charmonium spectrum seems 
academic. 

While the NRQCD approach appears to be problematic for the charmo- 
nium system, relativistic lattice formalisms have their share of difficulties. 
Simulations of charmonium with a variety of quark actions — NRQCD, the 
Fermilab actions, the D234 action — all underestimate the S- hyperfine split- 
ting by at least 40 MeV (see Q for a good summary). 

There are sound reasons for estimating the size of dynamical quark effects 
in the charmonium system. Some have suggested the remaining hyperfine 
discrepancy is due to quenching; estimates of the effects of dynamical quark 
loops range as high as 40% ||. Our results indicate this is unlikely to 
be the case — we find at most a 10% difference between our quenched and 
unquenched hyperfine splittings. As the quenching effects are apparently 
small for the range of different quark interactions present in the NRQCD 
action, we suggest that they will also be small across other quark actions. 

We recognise several shortcomings in our study: we have used different 
gluon actions for quenched and unquenched simulations, we have only ex- 
amined the effects of unquenching at a single dynamical quark mass and 
a single lattice spacing, and we have not attempted to extrapolate to the 
physical case of three light sea quark flavours. The first of these issues 
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was discussed in Section [O] above. To address the other objections with 
further simulations is beyond our present computational resources. In any 
case, such efforts are perhaps justified in simulations of the b quark, where 
systematic uncertainties are under better control and quenching effects are 
probably of comparable size to discretisation and radiative effects. For the 
charm system however, the much larger high-order relativistic errors and 
the large tadpole corrections dominate the effects of quenching. 

The sensitivity of the NRQCD corrections to the choice of tadpole factor 
is well established. This sensitivity should disappear with a higher-order 
treatment of the tadpole loops (and other radiative corrections) in lattice 
simulations. In practice, such a treatment is not yet available, and some 
choice for the tadpole factor is required. Our results add to the growing list 
of evidence if favour of calculating the tadpole correction factor from the 
mean link in the Landau gauge, in preference to the plaquette definition. 

The large effects we have encountered due to instabilities in the evolution 
equation should also be investigated further. These instabilities are doubt- 
less amplified in simulations of the charm quark, where the convergence of 
the NRQCD expansion is already questionable. Using an improved evolu- 
tion equation, as we have demonstrated, may bring the NRQCD approach 
into agreement with other quenched relativistic results for charmonium. 

Further, we have shown that 0(ot s ) radiative corrections may shift the 
spin-splittings by as much as 40%. While this is a crude estimate, the pos- 
sibility of such sizeable corrections in comparison with the small quenching 
effect gives us pause for consideration. Of particular note are unquenched 
results for the T spectrum in using the C(v 6 ) Hamiltonian, which in- 
dicate that remaining discrepancies with experiment are at the ten percent 
level — conceivably within the reach of radiative corrections. Perturbative 
calculations of the remaining radiative corrections to the NRQCD coeffi- 
cients, and those in other actions as well, will likely be necessary in the near 
future. 

We thank Howard Trottier and Randy Lewis for stimulating discussions. 
This work has been supported by the Natural Sciences and Engineering 
Research Council of Canada. 
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p 


v p 

u 




a (fm) a' 1 (GeV) 


aM 


M k 


s 








Quenched 








2.52 


0.874 




0.168(3) 1.17(2) 


0.81 


3.0(1) 


6 


2.10 




0.829 


0.181(3) 1.09(2) 


1.15 


3.0(1) 


4 








Unquenched 








5.415 


0.854 




0.163(3) 1.21(2) 


0.82 


2.9(1) 


6 


5.415 




0.800 


0.163(3) 1.21(2) 


1.15 


2.9(1) 


4 



Table 1: Parameters used in charmonium simulations. The lattice volume 
is 12 3 x 24 for the quenched simulations, and 16 3 x 32 for the unquenched 
simulations. is the kinetic mass of the 1 S state; s is the NRQCD 

stability parameter in Equation ([!]). 



tmin ■ t m ax 


l s 




3d I q 
Dl - Do 


2:24 


0.6646(4) 


0.7000(5) 




3:24 


0.6630(5) 


0.6984(5) 


0.03659(9) 


4:24 


0.6625(5) 


0.6979(6) 


0.03661(11) 


5:24 


0.6624(6) 


0.6977(7) 


0.03645(13) 


6:24 


0.6623(7) 


0.6976(7) 


0.0364(2) 


7:24 


0.6623(7) 


0.6976(8) 


0.0364(2) 


train ■ tmax 


'Pi 


3 Po 3 Pi 


3 P2 


2:14 


1.109(4) 


1.159(4) 1.138(5) 


1.082(4) 


3:14 


1.093(5) 


1.127(7) 1.113(6) 


1.072(5) 


4:14 


1.085(7) 


1.122(10) 1.102(9) 


1.065(6) 


5:14 


1.087(10) 


1.139(15) 1.107(12) 


1.067(9) 


6:14 


1.091(13) 


1.14(2) 1.114(18) 


1.067(12) 



Table 2: Examples of fits to quenched charmonium propagators. The fit 
results are shown for the plaquette-tadpole simulation, for various sets of 
{tmin '■ t m ax)- Single exponential fits are used for individual masses, and a 
ratio fit is used to extract the S-state hyperfine splitting. Italicised entries 
indicate the final results. 
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State 





0(v 4 ) 


0(v 6 ) 


0(v A ) 


0(v 6 ) 


1 Sq 


0.5733(5) 


0.6625(5) 


0.1708(4) 


0.2297(4) 


3 s 1 


0.6635(8) 


0.6979(6) 


0.2466(6) 


0.2802(5) 


'Pi 


1.034(8) 


1.085(7) 


0.643(7) 


0.696(7) 


3 P 


0.966(7) 


1.122(10) 


0.576(6) 


0.661(7) 


3 P1 


1.006(8) 


1.102(9) 


0.628(7) 


0.695(8) 


3 P2 


1.088(8) 


1.065(6) 


0.669(9) 


0.692(7) 


3 S 1 - 1 S 


0.0910(3) 


0.0365(1) 


0.0778(3) 


0.0521(2) 


State 


0(v 4 ) 


0{v 6 ) 


0{v A ) 0(v 6 ) 


Expt 



6 s, 

'Pi 



3 P 



3 P1 
3 P 



3.086(2) 
3.517(17) 
3.439(16) 
3.486(17) 
3.576(17) 
0.106(3) 



3.022(2) 
3.470(17) 
3.522(20) 
3.488(17) 
3.449(16) 
0.042(2) 



3.066(2) 
3.499(18) 
3.426(15) 
3.483(18) 
3.528(22) 
0.086(2) 



3.036(1) 
3.479(17) 
3.441(17) 
3.478(18) 
3.475(17) 
0.056(2) 



3.097 
3.526 
3.417 
3.511 
3.556 
0.118 



Table 3: Quenched charmonium masses in lattice units (top) and GeV 
(bottom). The scale is set by a -1 in Table |l[ 
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State 






«0 








0(v 6 ) 


0(v 4 ) 


0(v 6 ) 




0.5501(3) 


0.6279(3) 


0.0581(3) 


0.1155(3) 


3 s 1 


0.6363(5) 


0.6668(4) 


0.1278(4) 


0.1658(4) 


'Pi 


0.988(7) 


1.030(10) 


0.485(7) 


0.537(6) 


3 Po 


0.937(6) 


1.020(8) 


0.433(6) 


0.503(6) 


3 Pi 


0.977(7) 


1.050(10) 


0.476(7) 


0.531(7) 


3 P2 


1.016(9) 


1.065(5) 


0.497(9) 


0.540(7) 




0.0884(3) 


0.0365(2) 


0.0719(2) 


0.0522(1) 


State 






< 


Expt 




0(v 4 ) 


0{v 6 ) 


0{v A ) 0(v 6 ) 




3 Si 


3.087(2) 


3.028(2) 


3.068(2) 3.043(1) 


3.097 


l Px 


3.514(17) 


3.475(17) 


3.500(17) 3.486(16) 


3.526 


3 P 


3.456(15) 


3.518(20) 


3.437(15) 3.445(15) 


3.417 


3 Pi 


3.501(17) 


3.499(17) 


3.490(17) 3.479(17) 


3.511 


3 P2 


3.548(21) 


3.462(16) 


3.515(21) 3.489(17) 


3.556 


3 s 1 - 1 s 


0.107(2) 


0.049(1) 


0.087(2) 0.062(1) 


0.118 



Table 4: Unquenched charmonium masses in lattice units (top) and GeV 
(bottom). The scale is set by a -1 in Table g} 



5-state hyperfine P-state hyperfine 











3 P2~ 3 P0 


Experiment 






0.118 


0.139 


0{v A ) 


Quenched 


Uq 


0.106(2) 


0.14(2) 






^0 


0.086(2) 


0.10(2) 




Unquenched 




0.108(2) 


0.10(2) 






u% 


0.087(2) 


0.08(2) 


0(v 6 ) 


Quenched 


u 


0.042(1) 


-0.07(2) 






^0 


0.056(1) 


0.033(15) 




Unquenched 


u% 


0.049(1) 


-0.023(17) 






Uq 


0.062(1) 


0.044(15) 



Table 5: Quenched and unquenched 5-state and P-state hyperfine splittings 
for both Un and vk simulations. 
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Figure 1: The effective mass, — In ^ G q'(^ > ) , of the 1 Sq (circles) and 
(squares) in the quenched simulation. 



20 



1.4 



"I 1-2- 
=> 

O 

-I — ' 
-t— • 



w 



> 
o 

t 0.8 
LU 



0.6 



10 



12 



1 4 



Figure 2: The effective mass, — In ( ^q^)^ ' °^ ^ ne 3 ^ >0 m ^ ne quenched 
simulation. 
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Figure 3: Quenched charmonium spectrum using Uq. Squares represent 
results obtained to 0(v 4 ), circles represent 0(v 6 ) data. Horizontal lines 
indicate experimental values. 
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Figure 4: Quenched charmonium spectrum using Uq. Squares represent 
results obtained to 0(v ), circles represent 0(v 6 ) data. Horizontal lines 
indicate experimental values. 
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Figure 5: Unquenched charmonium spectrum using Uq . Squares represent 
results obtained to 0(v 4 ), circles represent 0(v 6 ) data. Horizontal lines 
indicate experimental values. 
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Figure 6: Unquenched charmonium spectrum using Uq . Squares represent 
results obtained to 0(v 4 ), circles represent 0(v 6 ) data. Horizontal lines 
indicate experimental values. 



25 



0.15 



O Quenched 
□ Unquenched 



> 
CD 



0.1 



Q. 

W 

CD 

CD 
Q- 



a) 0.05 

co 
-t— < 
w 

CO 



ffl 



Expt 



P4 



P6 



L4 



L6 



Figure 7: Charmonium 5-state hyperfine splitting. 'P4', 'P6' refer to the 
0(v 4 ), 0(v 6 ) results obtained with the plaquette tadpole factor; 'L4' and 
'L6' are the Landau tadpole results. 
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Figure 8: Charmonium P-state hyperfine splitting. 'P4', 'P6' refer to the 
0(v 4 ), 0(v 6 ) results obtained with the plaquette tadpole factor; 'L4' and 
'L6' are the Landau tadpole results. 
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